function oz_ngp_image,value,xdata_in,ydata_in,xlims_in,ylims_in,nx,ny,xl=xl,yl=yl

;+
; NAME:
;       OZ_NGP_IMAGE
;
; PURPOSE:
;       Assigns array values to a grid based on nearest grid point
;       scheme. Useful for things like T-rho planes with lots of particles
;
; CALLING SEQUENCE:
;       OZ_NGP_IMAGE(value,xpos,ypos,xlims,ylims,nx,ny,xl=xl,yl=yl)
;
; INPUTS:
;
; OPTIONAL INPUTS:
;
; OUTPUTS:
;
; MODIFICATION HISTORY:
;       Started by Rob Crain, Swinburne, 24-Feb-2010
;-

  ; declare image
  image = fltarr(nx,ny)

  ; convert to log if necessary 
  if keyword_set(xl) then begin
     xlims = alog10(xlims_in)
     xdata = alog10(xdata_in)
  endif else begin
     xlims = xlims_in
     xdata = xdata_in
  endelse

  if keyword_set(yl) then begin
     ylims = alog10(ylims_in)
     ydata = alog10(ydata_in)
  endif else begin
     ylims = ylims_in
     ydata = ydata_in
  endelse

  ; set step size
  dx = (xlims[1]-xlims[0]) / float(nx-1)
  dy = (ylims[1]-ylims[0]) / float(ny-1)

  ; assign grid indices to particles
  ix  = (xdata - xlims[0])/dx
  iy  = (ydata - ylims[0])/dy
  ind = where(ix GE 0 and ix LT nx and iy GE 0 and iy LT ny,n)
  if(n GT 0)then begin
     for i=long(0),long(n)-1l do image[ix[ind[i]],iy[ind[i]]] += value[ind[i]]
  endif

  return,image

end
